PEC3
Ferran Lugo y Yury Postnikov
12/22/2021
Ejercicio 1
- Buscar un conjunto de datos relacionados con la Bioestadística o la Bioinformática. Posibles recursos bibliográficos son aquellos presentados en la PEC1, por ejemplo: http://www.bioinformatics.org/sms2/index.html; o bien http://biostat.mc.vanderbilt.edu/wiki/Main/DataSets. También es posible utilizar fuentes propias o de interés, siempre teniendo en cuenta que sean datos públicos para los usuarios. Será necesario especificar la procedencia y la justificación de la elección de los datos.
Para nuestro trabajo hemos elegido la base de datos pública, que hemos encontrado en el sitio web kaggle.com:
https://www.kaggle.com/uciml/biomechanical-features-of-orthopedic-patients (data recodida inicialmente de Lichman, M. (2013). UCI Machine Learning Repository [http://archive.ics.uci.edu/ml]. Irvine, CA: University of California, School of Information and Computer Science)
Hemos eligido este data set por: - Una de las problemas frecuentes en la medicina en el ambito contemporaneo es la optimización de recursos donde sea posible. Y una área de aplicación de la analísis es en la etápa del diagnóstico de algunas enfermedades en las etapas tempranas, así por ejemplo, el diagnóstico de la hernia intervertebral exige utilizar la prubas caras y que tardan mucho en realizarse, cuando hay alternativas, que, quizas, no son tan fiables, pero puedan ayudar a la hora de selecionar pacientes para la siguiente etápa de estudios. En este caso se puede hablar de comparación entre estudio con tomografia computarizada vs la resonancia magnética. - Es un data-set público y facil de acceder. - Es un data-set relativamente límpio y ya preparado para realizar analisis de datos sobre el.
Ejercicio 2
- Utilizando R, mostrar y explicar qué tipo de ficheros se han importado, las variables de estudio (tipo, clasificación, …), así como todo aquello relevante para el estudio. Incluir capturas de pantalla y las instrucciones en R utilizadas para importar y mostrar los datos.
Esta dataset tiene 2 archivos .csv (column_2C_weka y column_3C_weka), que se puede juntar y valorar conjuntamente: Es una base de datos que contiene informacaión sobre 310 personas, en las que evaluaron posición de vertebras de la columna lumbar utilizando varios parametros (tabla 1) y exolicación visual con la imagen del articulo^1. Dentro de esta muestra habia 100 personas que consideramos “sanos” o “normales” y 210 (“anormales”) que consideramos que tenian algun tipo de patología: “hernia” o “spondilolistesis”.
tabla 1:
| Variable | Expicación del variable |
|---|---|
| pelvic_incidence, numerico en grados | Angulo entre la línea que pasa del centro de la lámina terminal del S1 al centro de las cabezas femorales, y la línea perpendicula a la lámina terminal de S1 |
| pelvic_tilt, numerico en grados | Angulo entre la línea que pasa del centro de la lámina terminal del S1 al centro de las cabezas femorales y una línea vertical |
| lumbar_lordosis_angle, numerico | angulo entre la lámina terminal superior del L1 y S1 |
| sacral_sclope, numerico en grados | angulo entre la línea perpendicula a la lámina terminal de S1 y una línea horizontal |
| pelvic_radius, numerico en grados | Es un angulo entre la línea vertical y una línea del centro de las cabezas femorales hasta el borde posterior de la lámina terminal de S1 |
| degree_spondylolysthesis, numerico en mm | distancia de desplazamiento de la vertebra superior sobre la vertebra inferior (positivo - anterior y negativo - posterior) |
| class | es posible utilizar este variable como character o factor. |
| class_2c, factor | Classificación de data como “normal” y “anoraml” |
| class_3c, factor | Classificación de data como “normal”, “Hernia” y “spondilolistesis” |
En esta imagen presentaremos la visualización de lo angulos medidas en la base de datos. 
Las bibliografía:
Sergides IG, McCombe PF, White G, Mokhtar S, Sears WR. Lumbo-pelvic lordosis and the pelvic radius technique in the assessment of spinal sagittal balance: strengths and caveats. Eur Spine J. 2011;20 Suppl 5(Suppl 5):591-601. doi:10.1007/s00586-011-1926-z
https://www.orthobullets.com/spine/2071/lumbar-spine-anatomy para la fecha de 25.12.2021
# Vamos a introducir los datos en RStudio. Para esto vamos a utilizar la libreria básica del R
column_2C <- read_csv("archive/column_2C_weka.csv",
col_types = cols(pelvic_incidence = col_number(),
`pelvic_tilt numeric` = col_number(),
lumbar_lordosis_angle = col_number(),
sacral_slope = col_number(), pelvic_radius = col_number(),
degree_spondylolisthesis = col_number()))
column_3C <- read_csv("archive/column_3C_weka.csv",
col_types = cols(pelvic_incidence = col_number(),
pelvic_tilt = col_number(), lumbar_lordosis_angle = col_number(),
sacral_slope = col_number(), pelvic_radius = col_number(),
degree_spondylolisthesis = col_number()))
# Vamos a ajustar los nombres de variables y confirmamos que son iguales
colnames(column_2C) <- colnames(column_3C)
colnames(column_2C) %>%
set_names() %>%
map(~(column_2C[,.x]==column_3C[,.x])) %>%
as.data.frame() %>%
summary() #utilizaremos función map de la libreria purr de tidyverse, aplicando la función a todos los elementos de las columnas correspondientes.## pelvic_incidence pelvic_tilt lumbar_lordosis_angle sacral_slope
## Mode:logical Mode:logical Mode:logical Mode:logical
## TRUE:310 TRUE:310 TRUE:310 TRUE:310
##
## pelvic_radius degree_spondylolisthesis class
## Mode:logical Mode:logical Mode :logical
## TRUE:310 TRUE:310 FALSE:210
## TRUE :100# Y cambiamos el nombre de la columna "class" para poder juntar los dos data frame.
column_2C<-rename(column_2C, "class_2" = "class")
column_3C<-rename(column_3C, "class_3" = "class")
# Creamos un data frame con información de los dos archivos iniciales.
ortho_df <- merge(column_2C, column_3C, by = c("pelvic_incidence", "pelvic_tilt", "lumbar_lordosis_angle", "sacral_slope", "pelvic_radius", "degree_spondylolisthesis"))
# Convertimos la varibel class en un factor, con los niveles que van de nivel normal a nivel con alteraciones más faciles de encontrar (no se puede evidenciar la hernía intervertebral en la radiografia, pero la espondilolistesis - sí).
ortho_df <- ortho_df %>%
mutate("class_2" = factor(ortho_df$`class_2`,
levels = c("Normal", "Abnormal")))
ortho_df <- ortho_df %>%
mutate("class_3" = factor(ortho_df$`class_3`,
levels = c("Normal", "Hernia", "Spondylolisthesis")))# En algunas ocasiones, para el objetivo de visualizar la data con mayor facilidad, vamos a cerar un data set con
ortho_df2 <- ortho_df %>%
pivot_longer(cols = c(pelvic_incidence, pelvic_tilt, pelvic_radius, sacral_slope, lumbar_lordosis_angle),
names_to = "type_of_measure", values_to = "angles") %>%
select(type_of_measure, angles, everything())
ortho_df3 <- ortho_df %>% filter(degree_spondylolisthesis<400) # limpiamos el data.frame de un valor evidentemente erroneo (degrre_spondilolistesis = 418 mm)Ejercicio 3:
- Con la Sección 2 de la PEC1 como base, elaborar y analizar una serie de cuestiones, que ayuden a explorar y a familiarizarse mejor con los datos de estudio. Además, en algunos casos, puede utilizarse la definición de funciones y el lenguaje SQL estudiado en el LAB3.
De todas las cuestiones que podemos realizar con estos datos, seleccionaremos las dos que, para nosotros, son las más destacables y extrapolables a una utilidad médica real.
La primera cuestión sería relativa a la comparación de grupos. Pregunta 1: ¿Existen diferencias significativas entre las variables angulares del grupo “Normal” respecto al grupo “Abnormal” (personas con hernia o con espondilolistesis)?. En otras palabras, ¿el hecho de tener hernia o espondilolistesis afecta significativamente en alguna de las variables de estudio? Su respuesta nos permitirá discernir qué métodos de medida difieren estadísticamente respecto al grupo “Normal” y llegar a un mejor acercamiento a nivel de diseño de prótesis o clasificación de características de pacientes en bases de datos.
Para resolver esta pregunta compararemos los datos de los grupos a lo largo de este estudio, visualizaremos los estadísticos básicos de nuestra muestra y realizaremos un análisis de varianza (ANOVA) de los grupos (véase los ejercicios 4 y 7).
Pregunta 2: ¿Existe alguna relación o correlación entre las variables del estudio? Esta pregunta nos permitirá conocer la relación entre las variables de estudio para saber si puede existir cierta redundancia en la medida de la anatomía de la columna lumbar. Además, nos permitirá saber si existe un patrón significativo de influencia de una variable determinada a otra/otras variable/s.
Esta cuestión podrá resolverse mediante los estudios de regresión líneal y regresión múltiple (véase los ejercicios 5 y 6). Además, se realizará un análisis de clústering (véase ejercicio 7) para observar si es posible agrupar los datos de nuestra muestra.
Ejercicio 4:
- Realizar un análisis descriptivo de los datos. El análisis debe incluir un resumen paramétrico de los datos y su representación gráfica, que mejor defina y complemente cada una de dichas variables.
El primer elemento de nuestros datos que analizaremos es la estadística básica y el comportamiento de cada variable entre los grupos.
Primeramente observaremos la estadística básica de cada variable mediante el comando “summary()”. Éste nos muestra los valores mínimos, máximos, primer y tercer cuartil, mediana y media para cada variable. Más adelante calcularemos su desviación estándar para cada grupo.
Al ser una variable de tipo caracter, mostramos las observaciones por grupos en las variables “class_2” y “class_3” mediante la herramienta “table()”.
summary(ortho_df)## pelvic_incidence pelvic_tilt lumbar_lordosis_angle sacral_slope
## Min. : 26.15 Min. :-6.555 Min. : 14.00 Min. : 13.37
## 1st Qu.: 46.43 1st Qu.:10.667 1st Qu.: 37.00 1st Qu.: 33.35
## Median : 58.69 Median :16.358 Median : 49.56 Median : 42.40
## Mean : 60.50 Mean :17.543 Mean : 51.93 Mean : 42.95
## 3rd Qu.: 72.88 3rd Qu.:22.120 3rd Qu.: 63.00 3rd Qu.: 52.70
## Max. :129.83 Max. :49.432 Max. :125.74 Max. :121.43
## pelvic_radius degree_spondylolisthesis class_2
## Min. : 70.08 Min. :-11.058 Normal :100
## 1st Qu.:110.71 1st Qu.: 1.604 Abnormal:210
## Median :118.27 Median : 11.768
## Mean :117.92 Mean : 26.297
## 3rd Qu.:125.47 3rd Qu.: 41.287
## Max. :163.07 Max. :418.543
## class_3
## Normal :100
## Hernia : 60
## Spondylolisthesis:150
##
##
## table(ortho_df$class_2)##
## Normal Abnormal
## 100 210table(ortho_df$class_3)##
## Normal Hernia Spondylolisthesis
## 100 60 150Obsérvese que podemos determinar los límites anotados para cada variable y el comportamiento general mediante la media. Estos valores nos servirán más adelante para determinar un intervalo para presentar las variables en gráficos.
La variable “pelvic_incidence” tiene valores de 26.15º hasta 129.83º.
La variable “pelvic_tilt” tiene valores de -6.555º hasta 49.432º.
La variable “lumbar_lordosis_angle” tiene valores de 14º hasta 125.74º.
La variable “sacral_slope” tiene valores de 13.37º hasta 121.43º.
La variable “pelvic_radius” tiene valores desde 70.08º hasta 163.07º.
La variable “degree_spondylolisthesis” tiene valores desde -11.058 mm hasta 418.543 mm.
Las variables “class_2 y class_3” incluyen 210 personas en el grupo “Abnormal” (60 con hernia y 150 con espondilolistesis) y 100 personas en el grupo “Normal”.
Aunque esta descripción nos sirve para acotar los valores de las variables con los que trabajamos, los datos comprenden todos grupos de estudio en conjunto. Nuestro objetivo principal es la comparación entre personas del grupo “Normal” (sin anomalías aparentes relativas a la postura pélvica y lumbar) y “Abnormal” (personas con hernia o espondilolistesis diagnosticada) para determinar si existen diferencias significativas en las variables.
El hecho de dividir los datos podría llegar a mostrar una diferencia en el comportamiento de los mismos. Es por eso que realizaremos un análisis de comparación de variables según los grupos. El valor que compararemos será la media de cada variable, con su desviación estándar anotada. Asimismo, para hacer un análisis visual, las diferencias entre grupos de estudio se representarán en gráficos de cajas (boxplots). En cada caso se utilizarán los intervalos calculados anteriormente para que la representación de cada gráfico esté en el mismo intervalo y que así su comparación sea más intuitiva y fácil.
Hace falta remarcar que este análisis lo haríamos para las variables de tipo angular (las cinco primeras variables). No obstante, al tratarse de un trabajo de ámbito educativo, para la variable “grado de espondilolistesis” realizaremos la comparación igualmente, pero de antemano sabemos que será superior en pacientes del grupo “espondilolistesis” respecto al resto, ya que tienen dicha afectación diagnosticada.
Mediante la herramienta “aggregate” calculamos la media y la desviación estándar para cada grupo de nuestras variables.
aggregate(ortho_df[, 1:6], list(ortho_df$class_3), mean)## Group.1 pelvic_incidence pelvic_tilt lumbar_lordosis_angle
## 1 Normal 51.68524 12.82141 43.54260
## 2 Hernia 47.63841 17.39879 35.46352
## 3 Spondylolisthesis 71.51422 20.74804 64.11011
## sacral_slope pelvic_radius degree_spondylolisthesis
## 1 38.86383 123.8908 2.186572
## 2 30.23961 116.4750 2.480251
## 3 50.76619 114.5188 51.896687aggregate(ortho_df[, 1:6], list(ortho_df$class_3), sd)## Group.1 pelvic_incidence pelvic_tilt lumbar_lordosis_angle
## 1 Normal 12.36816 6.778503 12.361388
## 2 Hernia 10.69713 7.016708 9.767795
## 3 Spondylolisthesis 15.10934 11.506169 16.397068
## sacral_slope pelvic_radius degree_spondylolisthesis
## 1 9.624004 9.014246 6.307483
## 2 7.555388 9.355720 5.531177
## 3 12.318813 15.579995 40.108030Representaremos estos estadísticos a continuación para cada variable.
Variable 1: Incidencia pélvica (pelvic_incidence)
par(mfrow=c(1,4))
boxplot(ortho_df$pelvic_incidence[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Incidencia pélvica", xlab="Normal", ylim=c(20,130))
boxplot(ortho_df$pelvic_incidence[ortho_df$class_2=="Abnormal"],col="Tomato",main="Incidencia pélvica", xlab="Abnormal",ylim=c(20,130))
boxplot(ortho_df$pelvic_incidence[ortho_df$class_3=="Hernia"],col="coral1",main="Incidencia pélvica", xlab="Hernia", ylim=c(20,130))
boxplot(ortho_df$pelvic_incidence[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Incidencia pélvica", xlab="Spondylolisthesis", ylim=c(20,130))
par(mfrow=c(1,1))Observamos que la media del grupo “Normal” es ligeramente inferior a la del grupo “Abnormal”. Al haber más pacientes en la muestra con “Spondylolisthesis” que en la de “Hernia”, la media del grupo “Abnormal” resulta ser superior a la del grupo “Normal”, con una desviación estándar muy alta. Este hecho se observa en los gráficos, donde parece que el grupo “Spondylolisthesis” supera significativamente la media otorgada por el grupo “Normal”. Como sólo observamos tres outliers aparentes en el grupo de espondilolistesis (la muestra no parece tener mucha distorsión), nuestra primera hipótesis será que los pacientes con espondilolistesis podrían tener una mayor incidencia pélvica general. Así, es probable que la variable “incidencia pélvica” tenga cierto peso a la hora de analizarse en dicho grupo.
Variable 2: Inclinación pélvica (pelvic_tilt)
par(mfrow=c(1,4))
boxplot(ortho_df$pelvic_tilt[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Inclinación pélvica", xlab="Normal", ylim=c(-7,50))
boxplot(ortho_df$pelvic_tilt[ortho_df$class_2=="Abnormal"],col="Tomato",main="Inclinación pélvica", xlab="Abnormal",ylim=c(-7,50))
boxplot(ortho_df$pelvic_tilt[ortho_df$class_3=="Hernia"],col="coral1",main="Inclinación pélvica", xlab="Hernia", ylim=c(-7,50))
boxplot(ortho_df$pelvic_tilt[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Inclinación pélvica", xlab="Spondylolisthesis", ylim=c(-7,50))
par(mfrow=c(1,1))Se observa que la media del grupo “Normal” para la inclinación pélvica es ligeramente inferior que la de los otros grupos. No obstante, al existir una desviación estándar realmente alta para todos los grupos en esta variable (es decir, existe una gran variabilidad), es probable que la distinción entre grupos no llegue a ser significativa. Por lo tanto, no esperamos diferencias significativas entre individuos de los distintos grupos para esta variable.
Variable 3: Ángulo de lordosis lumbar (lumbar_lordosis_angle)
par(mfrow=c(1,4))
boxplot(ortho_df$lumbar_lordosis_angle[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Lordosis lumbar", xlab="Normal", ylim=c(13,126))
boxplot(ortho_df$lumbar_lordosis_angle[ortho_df$class_2=="Abnormal"],col="Tomato",main="Lordosis lumbar", xlab="Abnormal",ylim=c(13,126))
boxplot(ortho_df$lumbar_lordosis_angle[ortho_df$class_3=="Hernia"],col="coral1",main="Lordosis lumbar", xlab="Hernia", ylim=c(13,126))
boxplot(ortho_df$lumbar_lordosis_angle[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Lordosis lumbar", xlab="Spondylolisthesis", ylim=c(13,126))
par(mfrow=c(1,1))En este caso observamos que la media del grupo “Normal” es parecida a la del grupo “Hernia”, pero difiere bastante de la del grupo “Spondylolisthesis”. Aunque exista algún outlier apartado que pueda hacer aumentar la desviación estándar de cada grupo (véase los puntos superiores que se distancian de la media en cada gráfico), es probable que esta comparativa pueda ser significativa. Así, nuestra segunda hipótesis es que los pacientes con espondilolistesis podrían tener un mayor ángulo de lordosis lumbar. Y por lo tanto, la variable “ángulo de lordosis lumbar” podría tener más peso en este grupo de pacientes.
Variable 4: Pendiente sacral (sacral_slope)
par(mfrow=c(1,4))
boxplot(ortho_df$sacral_slope[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Pendiente sacral", xlab="Normal", ylim=c(12,122))
boxplot(ortho_df$sacral_slope[ortho_df$class_2=="Abnormal"],col="Tomato",main="Pendiente sacral", xlab="Abnormal",ylim=c(12,122))
boxplot(ortho_df$sacral_slope[ortho_df$class_3=="Hernia"],col="coral1",main="Pendiente sacral", xlab="Hernia", ylim=c(12,122))
boxplot(ortho_df$sacral_slope[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Pendiente sacral", xlab="Spondylolisthesis", ylim=c(12,122))
par(mfrow=c(1,1))Para la pendiente sacral observamos que las medias de los grupos difieren muy ligeramente, siendo superior en el grupo “Spondylolisthesis”. Existen outliers que pueden alterar la media de los datos, sobretodo en el grupo “Spondylolisthesis”, donde se observan varios puntos que se alejan significativamente del comportamiento del resto de las muestras. Es por eso que probablemente no existan diferencias significativas entre la media de los grupos para esta variable, ya que sin estos outliers la media del grupo “Spondylolisthesis” podría disminuir y se asemejaría aún más a la de los otros grupos.
Variable 5: Radio pélvico (pelvic_radius)
par(mfrow=c(1,4))
boxplot(ortho_df$pelvic_radius[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Radio pélvico", xlab="Normal", ylim=c(69,164))
boxplot(ortho_df$pelvic_radius[ortho_df$class_2=="Abnormal"],col="Tomato",main="Radio pélvico", xlab="Abnormal",ylim=c(69,164))
boxplot(ortho_df$pelvic_radius[ortho_df$class_3=="Hernia"],col="coral1",main="Radio pélvico", xlab="Hernia", ylim=c(69,164))
boxplot(ortho_df$pelvic_radius[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Radio pélvico", xlab="Spondylolisthesis", ylim=c(69,164))
par(mfrow=c(1,1))Las medias de todos los grupos son semejantes, siendo ligeramente inferiores aquellas del grupo abnormal (incluyendo a los grupos hernia y espondilolistesis por separado). La gran variabilidad en el grupo de espondilolistesis (fíjese en su gran desviación estándar) hace dudar de la significación de esta comparación, aunque sea posible que la variable “radio pélvico” sea distintiva en el grupo “Normal” respecto al resto. Aún así, debido a la variabilidad no esperamos que se cumpla esta hipótesis a priori.
Variable 6: Grado de espondilolistesis (degree_spondylolisthesis)
par(mfrow=c(1,4))
boxplot(ortho_df$degree_spondylolisthesis[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="G.espondilolistesis", xlab="Normal", ylim=c(-12,419))
boxplot(ortho_df$degree_spondylolisthesis[ortho_df$class_2=="Abnormal"],col="Tomato",main="G.espondilolistesis", xlab="Abnormal",ylim=c(-12,419))
boxplot(ortho_df$degree_spondylolisthesis[ortho_df$class_3=="Hernia"],col="coral1",main="G.espondilolistesis", xlab="Hernia", ylim=c(-12,419))
boxplot(ortho_df$degree_spondylolisthesis[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="G.espondilolistesis", xlab="Spondylolisthesis", ylim=c(-12,419))
par(mfrow=c(1,1))Efectivamente el grado de espondilolistesis en pacientes con espondilolistesis es superior respecto al resto, el cual varía ligeramente de la media 0. Esta comparación hace clara la presencia de un elemento de la muestra probablemente mal tomado (existe un paciente con 400 mm (40 cm) de desplazamiento vertebral, hecho bastante cuestionable). No obstante, la comparativa sigue vigente aun ignorar los outliers que distorsionan la media.
Por otro lado, fíjese que esta variable es la única que se mide en milímetros, mostrándonos que la media de variación de la distancia entre vértebras lumbares en pacientes con espondilolistesis es de 51.89 mm. Por lo tanto, la mayoría de pacientes con espondilolistesis del estudio tienen un grado de espondilolistesis positivo, es decir, la vértebra desplazada tiene un desplazamiento anterior respecto su vértebra inferior en la mayoría de casos.
Ejercicio 5:
- Complementando el apartado anterior, elaborar un análisis de regresión de dos conjuntos de variables (LAB2 y Ejercicio 6 de la PEC1). La elección de las variables, los resultados, así como su relación deben de estar correctamente justificada.
Para realizar la regresión lineal eligiremos las dos variables que parezcan más relacionadas de entre todas las disponibles. Para ver el comportamiento que cada variable tiene con las otras variables del estudio, utilizaremos un gráfico de dispersión de pares de variables (pairs). Nótese que eliminaremos del análisis las dos variables de tipo caracter (class_2 y class_3) debido a que no son numéricas. Para analizar la relación entre variables de forma cuantitativa, calcularemos las correlaciones entre las variables numéricas mediante el comando “cor()”.
pairs(ortho_df[, -c(7:8)])
cor(ortho_df[, -c(7:8)])## pelvic_incidence pelvic_tilt lumbar_lordosis_angle
## pelvic_incidence 1.0000000 0.62919877 0.71728236
## pelvic_tilt 0.6291988 1.00000000 0.43276386
## lumbar_lordosis_angle 0.7172824 0.43276386 1.00000000
## sacral_slope 0.8149600 0.06234529 0.59838689
## pelvic_radius -0.2474672 0.03266781 -0.08034361
## degree_spondylolisthesis 0.6387427 0.39786228 0.53366701
## sacral_slope pelvic_radius degree_spondylolisthesis
## pelvic_incidence 0.81495999 -0.24746721 0.63874275
## pelvic_tilt 0.06234529 0.03266781 0.39786228
## lumbar_lordosis_angle 0.59838689 -0.08034361 0.53366701
## sacral_slope 1.00000000 -0.34212835 0.52355746
## pelvic_radius -0.34212835 1.00000000 -0.02606501
## degree_spondylolisthesis 0.52355746 -0.02606501 1.00000000# Para que sea más facil interpretar las relaciones entre las variables vamos a crear una visualización de las interacciones en contexto de classificación Normal-Anormal y Normal-Hernia-Espondilolistesis.
model.matrix(~0. + pelvic_incidence+pelvic_tilt+pelvic_radius+sacral_slope+lumbar_lordosis_angle+ degree_spondylolisthesis + class_3, ortho_df) %>%
cor() %>%
ggcorrplot(lab=T,
lab_size=2.5,
type = "lower")
model.matrix(~0. + pelvic_incidence+pelvic_tilt+pelvic_radius+sacral_slope+lumbar_lordosis_angle+ degree_spondylolisthesis + class_3, ortho_df) %>%
cor() %>%
ggcorrplot(lab=T,
lab_size=2.5,
type = "lower")
Este comando (cor()) nos muestra los coeficientes de correlación entre cada variable, por defecto mediante el método de Pearson. Cuando más próximos a 1 o -1, más relación tienen las dos variables, y por lo tanto, el incremento (en caso de coeficiente positivo) o el decremento (en caso de coeficiente negativo) de una variable influencia al incremento/decremento de la otra variable, respectivamente.
Elección de variables De todos los coeficientes de correlación de Pearson entre dos variables distintas, los que tienen la cifra más elevada en valor absoluto son la incidencia pélvica (pelvic_incidence) con la pendiente sacral (sacral_slope), con un coeficiente de 0.8149600; seguido de la incidencia pélvica (pelvic_incidence) con el ángulo de lordosis lumbar (lumbar_lordosis_angle), con con un coeficiente de correlación de 0.71728236.
Hace falta remarcar que existen algunas asociaciones de variables más con coeficientes de correlación elevados que quizás sean importantes en el siguiente análisis de regresión múltiple (véase el Ejercicio 6). No obstante, en este ejercicio, queremos comprobar si es cierto que la incidencia pélvica es una variable que se correlaciona con la pendiente sacral. Al tener un coeficiente de correlación tan próximo a 1, nuestra hipósis es que sí que se cumple esta relación.
De las comparaciones visuales anteriores, si nos fijamos en el gráfico de comparación entre estas dos variables, observamos cómo la nube de puntos de los datos de nuestra muestra forma una línea en diagonal ascendente. Este hecho verifica que el coeficiente de correlación sea tan elevado y que éste tenga un valor positivo: cuando una de las variables aumenta, parece ser que la otra también.
Modelo de regresión lineal Para determinar si esta relación es significativa realizaremos un análisis de regresión lineal entre estas dos variables. Para ello utilizaremos un modelo de regresión lineal (lm) que denominaremos como “ortho_lm”. Nuestra hipótesis será que la pendiente sacrial es la variable independiente (explicativa) y que la incidencia pélvica es la variable dependiente (explicada), ya que el hecho de tener una pendiente sacral más elevada implicaría que el sacro de la persona en cuestión se inclinara hacia atrás y por lo tanto ésta tendría posiblemente un ángulo de incidencia pélvica superior.
ortho_lm <-lm(ortho_df$pelvic_incidence~ortho_df$sacral_slope)
summary(ortho_lm)##
## Call:
## lm(formula = ortho_df$pelvic_incidence ~ ortho_df$sacral_slope)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.359 -6.566 -1.403 4.383 32.989
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.5461 1.9079 8.148 9.35e-15 ***
## ortho_df$sacral_slope 1.0465 0.0424 24.680 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.01 on 308 degrees of freedom
## Multiple R-squared: 0.6642, Adjusted R-squared: 0.6631
## F-statistic: 609.1 on 1 and 308 DF, p-value: < 2.2e-16Mediante el “summary” del modelo observamos distintos valores estadísticos. El valor principal que nos interesa es el coeficiente de determinación (el R-squared) del modelo, que nos dictará qué tan cerca están los datos de una regresión lineal. En este caso, el coeficiente de determinación no es muy próximo a 1 o -1 (es 0.6642) por lo que la bondad del modelo no parece tan fuerte como creíamos. Según este coeficiente, este modelo explicaría tan solo un 66.42% de la variabilidad de los datos.
Normalidad y varianza de los errores Para comprobar si el coeficiente obtenido es válido, analizaremos la normalidad de los errores del modelo mediante un histograma básico, un boxplot y el gráfico de cuartiles de los residuos.
residuos <-rstandard(ortho_lm)
par(mfrow=c(1,3))
hist(residuos, main = "Histograma de residuos")
boxplot(residuos, main = "Diagrama de cajas de residuos")
qqnorm(residuos, main = "G. cuartiles de residuos")
qqline(residuos)
par(mfrow=c(1,1))En el primer gráfico observamos cómo, a nivel visual, los residuos parecen tener un comportamiento normal (la forma de una pendiente simétrica con un pico próximo al cero). El segundo gráfico nos añade la información relativa a los outliers. Parece ser que existen diversos valores de nuestra muestra que hacen alejarse los errores de la normalidad de los datos. No obstante, según el tercer gráfico, no parece ser significativo, ya que los residuos siguen una recta diagonal con una desviación en el tercer cuartil únicamente. Así, podemos concluir que los residuos siguen una distribución normal.
Una vez verificada la normalidad de los residuos, analizaremos si la varianza de los errores es constante. El hecho de que la varianza de los errores sea constante implicaría que la variable dependiente variaría debido a nuestra variable independiente, no debido a otros factores externos a nuestro análisis. Para ello, compararemos el gráfico de nuestro modelo con los residuos estandarizados.
plot(fitted.values(ortho_lm),rstandard(ortho_lm),xlab="Valores ajustados",ylab="Residuos estandarizados", col="darkorchid3")
abline(h=0)
Verificamos que existen algunas anotaciones de nuestra muestra (outliers) que se alejan del comportamiento de la mayoría (es decir, que se alejan de la línea centrada al cero). Aun así, éstas parecen no ser significativas debido a que son pocas.
Es por eso que aceptamos el modelo de regresión lineal y suponemos que el hecho de que su coeficiente de determinación no llegue a ser tan elevado como su coeficiente de correlación podría ser debido a la influencia de estos outliers.
Modelo final y conclusiones Finalmente, creemos que el modelo estudiado parece válido y lo representaremos en el siguiente gráfico. Véase cómo la mayoría de datos de nuestra muestra para las dos variables sigue un mismo comportamiento lineal ascendente.
plot(ortho_df$sacral_slope,ortho_df$pelvic_incidence, main="Modelo de regresión lineal", xlab="Pendiente sacral", ylab= "Incidencia pélvica")
abline(ortho_lm)
Como conclusión podemos decir que, aunque el coeficiente de determinación no sea muy elevado, las dos variables analizadas (incidencia pélvica y pendiente sacral) tienen cierta relación en cuanto a su comportamiento. Además, sabemos que esta relación es positiva (un aumento de la pendiente sacral puede implicar un aumento de la incidencia pélvica). Estos datos son consistentes con el comportamiento biológico esperado, ya que el hecho de tener una pendiente sacral más elevada significa que el ángulo entre el sacro y el eje horizontal aumenta, haciendo inclinar así la pelvis. Al hacerlo, el ángulo de incidencia pélvica aumentaría debido a que éste tiene su base de medida entre el eje bicoxofemoral (que no habría cambiado) y línea perpendicular a la placa sacra (que se habría movido).
Este argumento se sustenta de bibliografía, como el estudio de E.V. Geiger et al. (2007), que econtró una correlación entre estas dos variables en su muestra de pacientes.
Ejercicio 6:
- Realizar, a partir de los conceptos trabajados en el LAB4 y la PEC2, un estudio probabilístico (a elección propia) de al menos 3 de las variables, que ayude a esclarecer cuestiones de relevancia que se pl antean en los ámbitos de acción estudiados.
Con el análisis del ejercicio 5 hemos podido responder parcialmente la pregunta 2 que nos proponíamos: ¿existe alguna relación entre nuestras variables de estudio? Mediante la regresión lineal hemos comprobado que existe una correlación positiva significativa entre la pendiente sacral y la incidencia pélvica. No obstante, también hemos visto unos índices de correlación elevados entre otras variables de estudio. Es por eso que ahora realizaremos una regresión múltiple para analizar qué variables tienen relación directa con otra/s.
Elección de variables Para definir un modelo inicial de regresión múltiple, debemos indagar de nuevo en los coeficientes de correlación de nuestras variables, para así elegir aquellos que sean próximos a 1 o -1. Nótese que trabajaremos con el mismo método de Pearson para calcularlos.
cor(ortho_df[, -c(7:8)])## pelvic_incidence pelvic_tilt lumbar_lordosis_angle
## pelvic_incidence 1.0000000 0.62919877 0.71728236
## pelvic_tilt 0.6291988 1.00000000 0.43276386
## lumbar_lordosis_angle 0.7172824 0.43276386 1.00000000
## sacral_slope 0.8149600 0.06234529 0.59838689
## pelvic_radius -0.2474672 0.03266781 -0.08034361
## degree_spondylolisthesis 0.6387427 0.39786228 0.53366701
## sacral_slope pelvic_radius degree_spondylolisthesis
## pelvic_incidence 0.81495999 -0.24746721 0.63874275
## pelvic_tilt 0.06234529 0.03266781 0.39786228
## lumbar_lordosis_angle 0.59838689 -0.08034361 0.53366701
## sacral_slope 1.00000000 -0.34212835 0.52355746
## pelvic_radius -0.34212835 1.00000000 -0.02606501
## degree_spondylolisthesis 0.52355746 -0.02606501 1.00000000En este caso observamos que, a parte de la correlación significativa entre la incidencia pélvica y la pendiente sacral observada anteriormente, existen asociaciones con coeficientes de correlación elevados.
Es el caso de la incidencia pélvica (pelvic_incidence) con el ángulo de lordosis lumbar (lumbar_lordosis_angle) que tienen un coeficiente de correlación de 0.71728236. Seguidamente, tenemos la incidencia pélvica con el grado de espondilolistesis (degree_spondylolisthesis) (coeficiente de 0.6387427) y la incidencia pélvica con la inclinación pélvica (pelvic_tilt) (coeficiente de 0.6291988). Finalmente existe un coeficiente de correlación de 0.59838689 entre el ángulo de lordosis lumbar y la pendiente sacral (sacral_slope), un coeficiente de 0.53366701 entre el ángulo de lordosis lumbar y el grado de espondilolistesis, y un coeficiente de 0.52355746 entre la pendiente sacral y el grado de espondilolistesis.
A partir de estos casos, tenemos coeficientes de correlación en valor absoluto inferiores a 0.5 por lo que no los incluiremos a priori como destacables.
Por lo tanto, observamos que la variable de incidencia pélvica es la que parece tener más relación con las otras, teniendo los coeficientes de correlación de Pearson más elevados para cuatro de las otras variables. Es por eso que parece intuitivo que esta variable será nuestra variable objetivo para comparar con las otras. En este caso, incluiremos todas las otras variables en el modelo, aunque sabemos que probablemente el radio pélvico sea descartado más adelante, ya que es la única variable de nuestro estudio sin un coeficiente de correlación elevado.
Modelo de regresión múltiple Realizamos el modelo de regresión de la incidencia pélvica con las cinco variables numéricas restantes.
ortho_modelomu<-lm(pelvic_incidence~pelvic_tilt+lumbar_lordosis_angle+sacral_slope+pelvic_radius+degree_spondylolisthesis, data=ortho_df)
summary(ortho_modelomu)##
## Call:
## lm(formula = pelvic_incidence ~ pelvic_tilt + lumbar_lordosis_angle +
## sacral_slope + pelvic_radius + degree_spondylolisthesis,
## data = ortho_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.065e-08 -4.684e-10 -9.070e-11 4.062e-10 1.100e-08
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.423e-10 3.168e-09 -1.400e-01 0.889
## pelvic_tilt 1.000e+00 3.263e-11 3.065e+10 <2e-16 ***
## lumbar_lordosis_angle 1.887e-11 2.104e-11 8.970e-01 0.371
## sacral_slope 1.000e+00 3.038e-11 3.292e+10 <2e-16 ***
## pelvic_radius 6.219e-12 2.192e-11 2.840e-01 0.777
## degree_spondylolisthesis -5.313e-12 9.453e-12 -5.620e-01 0.575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.688e-09 on 304 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 8.354e+20 on 5 and 304 DF, p-value: < 2.2e-16El “summary” del modelo nos informa de la significación de la correlación de nuestra variable objetivo (incidencia pélvica) con cada otra variable. Obsérvese que los coeficientes significativos son únicamente con la inclinación pélvica (pelvic_tilt), con p-valor inferior a 2e-16, y con la pendiente sacral (sacral_slope), con p-valor inferior a 2e-16. Los asteriscos nos informan visualmente de que éstas son las únicas dos variables con correlaciones significativas para la incidencia pélvica.
Además, nos otorga el coeficiente de determinación (R-squared) del modelo, siendo éste el máximo (1). Esta información parece sospechosa, pues en biología es difícil obtener un coeficiente de determinación total, pero si nos fijamos en el p-valor del modelo, observamos que éste es significativo (inferior a 2.2e-16). Es probable que tengamos un caso de multicolinealidad, es decir, varias variables que tienen correlación muy elevada, y eso haría incrementar el coeficiente de determinación. Aun así, como es un análisis significativo, ésto nos verifica que la bondad del modelo de regresión es positiva. En otras palabras, toda la variabilidad de la variable en cuestión (incidencia pélvica) viene dada por otras variables del estudio.
Ahora que sabemos que nuestras variables predictoras para la incidencia pélvica son la pendiente sacral y la inclinación pélvica, perfeccionaremos el modelo de regresión.
ortho_modelomu2<-lm(pelvic_incidence~pelvic_tilt+sacral_slope, data=ortho_df)
summary(ortho_modelomu2)##
## Call:
## lm(formula = pelvic_incidence ~ pelvic_tilt + sacral_slope, data = ortho_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.038e-08 -3.259e-10 -7.410e-11 4.869e-10 1.077e-08
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.766e-10 9.827e-10 7.900e-01 0.43
## pelvic_tilt 1.000e+00 2.662e-11 3.757e+10 <2e-16 ***
## sacral_slope 1.000e+00 1.985e-11 5.039e+10 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.674e-09 on 307 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 2.101e+21 on 2 and 307 DF, p-value: < 2.2e-16En este sentido, si aceptamos el modelo, la variable “incidencia pélvica” tendría toda su variabilidad explicada por la pendiente sacral (como hemos corroborado en el ejercicio anterior) y la inclinación pélvica.
A nivel anatómico, estos datos tienen mucho sentido. Como se ha visto en el ejercicio anterior, existe una relación documentada entre la incidencia pélvica y la pendiente sacral. En cuanto a la incidencia pélvica y la inclinación pélvica, pasa exactamente lo mismo.
De hecho, tal como indica el artículo de J.C.Le Huec et al. (2011), en medicina la incidencia pélvica (PI) se calcula mediante la suma de la pendiente sacral (SS) más la inclinación pélvica (PT).
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3175921/
Por lo tanto, parece lógico que estas tres variables tengan una correlación positiva significativa y que, por lo tanto, a mayor inclinación pélvica, mayor sea la incidencia pélvica.
Ejerecicio 7:
- Complementando el apartado anterior, elaborar un análisis ANOVA de dos conjuntos de variables (LAB5 y Ejercicio 6 de la PEC2). La elección de las variables, los resultados, así como su relación deben de estar correctamente justificada. Además, realizar una prueba cluster de las variables, y si existe un fuerte agrupamiento, elaborar un dendograma (LAB5).
Para decidir si se puede utuilizar el test ANOVA, tenemos que confirmar inialmente si las variables sigue la distribución normal y si la variable no es homoscedastica (tiene la misma varianza).
# Empezamos mmirando a las variables y calculando valores medias y error estandart.
ortho_df2 %>% group_by(type_of_measure) %>% summarise("Mean of angle" = mean(angles), "Standart deviation" = sd(angles))## # A tibble: 5 × 3
## type_of_measure `Mean of angle` `Standart deviation`
## <chr> <dbl> <dbl>
## 1 lumbar_lordosis_angle 51.9 18.6
## 2 pelvic_incidence 60.5 17.2
## 3 pelvic_radius 118. 13.3
## 4 pelvic_tilt 17.5 10.0
## 5 sacral_slope 43.0 13.4ortho_df2 %>% summarise("Mean of displacement" = mean(degree_spondylolisthesis), "SD of displacement" = sd(degree_spondylolisthesis))## # A tibble: 1 × 2
## `Mean of displacement` `SD of displacement`
## <dbl> <dbl>
## 1 26.3 37.5# Creamos la visualización de la distribución de las variables con denisity plot y con boxplot, marcando valores médias. En la ditribución normal valor media y valor mediana y moda coinciden. Por la visualización, parece que la variable pelvic_radius y la sacral_slope pueden seguir la distribución normal. También miramos si el hecho de pertenecer a una otra classificación puede afectar drasticamente al tipo de la sistribución. Se evidencia el cambio importante de comportamiento de la variable degree_spondilolistesis en comparación normal/anormal y espondilolistesis vs el resto de classificaciones.
pt <- ortho_df2 %>% ggplot(aes(angles, fill = type_of_measure)) + geom_density(alpha = 0.2)
ptc2 <- ortho_df2 %>% ggplot(aes(angles, fill = type_of_measure)) + geom_density(alpha = 0.2)+ facet_grid(~class_2)
ptc3 <- ortho_df2 %>% ggplot(aes(angles, fill = type_of_measure)) + geom_density(alpha = 0.2)+ facet_grid(~class_3)
pp <- list(pt,ptc2,ptc3)
ggarrange(nrow = 3, ncol= 1, plotlist = pp)
pds <- ortho_df2 %>% ggplot(aes(degree_spondylolisthesis)) + geom_density(alpha = 0.2, fill = "darkgreen")
pds2 <- ortho_df2 %>% ggplot(aes(degree_spondylolisthesis, fill = class_2)) + geom_density(alpha = 0.2)+ facet_grid(~class_2)
pds3 <- ortho_df2 %>% ggplot(aes(degree_spondylolisthesis, fill = class_3)) + geom_density(alpha = 0.2)+ facet_grid(~class_3)
ppds <- list(pds,pds2,pds3)
ggarrange(nrow = 3, ncol= 1, plotlist = ppds)
# Evaluamos las valores medias y del error estandar segun la classificación:
# Segun class_2
ortho_df2 %>% group_by(type_of_measure, class_2) %>%
summarise("Mean of angle" = mean(angles), "Standart deviation" = sd(angles))## `summarise()` has grouped output by 'type_of_measure'. You can override using the `.groups` argument.## # A tibble: 10 × 4
## # Groups: type_of_measure [5]
## type_of_measure class_2 `Mean of angle` `Standart deviation`
## <chr> <fct> <dbl> <dbl>
## 1 lumbar_lordosis_angle Normal 43.5 12.4
## 2 lumbar_lordosis_angle Abnormal 55.9 19.7
## 3 pelvic_incidence Normal 51.7 12.4
## 4 pelvic_incidence Abnormal 64.7 17.7
## 5 pelvic_radius Normal 124. 9.01
## 6 pelvic_radius Abnormal 115. 14.1
## 7 pelvic_tilt Normal 12.8 6.78
## 8 pelvic_tilt Abnormal 19.8 10.5
## 9 sacral_slope Normal 38.9 9.62
## 10 sacral_slope Abnormal 44.9 14.5ortho_df2 %>% group_by(class_2) %>%
summarise("Mean of displacement" = mean(degree_spondylolisthesis), "SD of displacement" = sd(degree_spondylolisthesis))## # A tibble: 2 × 3
## class_2 `Mean of displacement` `SD of displacement`
## <fct> <dbl> <dbl>
## 1 Normal 2.19 6.28
## 2 Abnormal 37.8 40.6# Según class_3
ortho_df2 %>% group_by(type_of_measure, class_3) %>%
summarise("Mean of angle" = mean(angles), "Standart deviation" = sd(angles))## `summarise()` has grouped output by 'type_of_measure'. You can override using the `.groups` argument.## # A tibble: 15 × 4
## # Groups: type_of_measure [5]
## type_of_measure class_3 `Mean of angle` `Standart deviation`
## <chr> <fct> <dbl> <dbl>
## 1 lumbar_lordosis_angle Normal 43.5 12.4
## 2 lumbar_lordosis_angle Hernia 35.5 9.77
## 3 lumbar_lordosis_angle Spondylolisthesis 64.1 16.4
## 4 pelvic_incidence Normal 51.7 12.4
## 5 pelvic_incidence Hernia 47.6 10.7
## 6 pelvic_incidence Spondylolisthesis 71.5 15.1
## 7 pelvic_radius Normal 124. 9.01
## 8 pelvic_radius Hernia 116. 9.36
## 9 pelvic_radius Spondylolisthesis 115. 15.6
## 10 pelvic_tilt Normal 12.8 6.78
## 11 pelvic_tilt Hernia 17.4 7.02
## 12 pelvic_tilt Spondylolisthesis 20.7 11.5
## 13 sacral_slope Normal 38.9 9.62
## 14 sacral_slope Hernia 30.2 7.56
## 15 sacral_slope Spondylolisthesis 50.8 12.3ortho_df2 %>% group_by(class_3) %>%
summarise("Mean of displacement" = mean(degree_spondylolisthesis), "SD of displacement" = sd(degree_spondylolisthesis))## # A tibble: 3 × 3
## class_3 `Mean of displacement` `SD of displacement`
## <fct> <dbl> <dbl>
## 1 Normal 2.19 6.28
## 2 Hernia 2.48 5.49
## 3 Spondylolisthesis 51.9 40.0# Creamos la visualización de valores médias y error estandar de las variables.
angles_ms <- ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure, color = type_of_measure)) +
stat_summary(fun = mean,
shape = 18,
size = 1)+
stat_summary(fun = mean,
fun.min = function(x) {mean(x) - sd(x)},
fun.max =function(x) {mean(x) + sd(x)},
geom = "errorbar")+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())
angles_ms2 <- ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure, color = type_of_measure)) +
stat_summary(fun = mean,
shape = 18,
size = 1)+
stat_summary(fun = mean,
fun.min = function(x) {mean(x) - sd(x)},
fun.max =function(x) {mean(x) + sd(x)},
geom = "errorbar")+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank()) +
facet_wrap(~class_2)
angles_ms3 <- ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure, color = type_of_measure)) +
stat_summary(fun = mean,
shape = 18,
size = 1)+
stat_summary(fun = mean,
fun.min = function(x) {mean(x) - sd(x)},
fun.max =function(x) {mean(x) + sd(x)},
geom = "errorbar")+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank()) +
facet_wrap(~class_3)
degree_ms <- ortho_df2 %>% ggplot(aes(y= degree_spondylolisthesis, x = 0)) +
stat_summary(fun = mean,
shape = 18,
size = 1,
color = "green")+
stat_summary(fun = mean,
fun.min = function(x) {mean(x) - sd(x)},
fun.max =function(x) {mean(x) + sd(x)},
color = "green",
geom = "errorbar")+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())
degree_ms2 <- ortho_df2 %>% ggplot(aes(y= degree_spondylolisthesis, x = class_2, color = class_2)) +
stat_summary(fun = mean,
shape = 18,
size = 1)+
stat_summary(fun = mean,
fun.min = function(x) {mean(x) - sd(x)},
fun.max =function(x) {mean(x) + sd(x)},
geom = "errorbar")+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank()) +
facet_wrap(~class_2)
degree_ms3 <- ortho_df2 %>% ggplot(aes(y= degree_spondylolisthesis, x = class_3, color = class_3)) +
stat_summary(fun = mean,
shape = 18,
size = 1)+
stat_summary(fun = mean,
fun.min = function(x) {mean(x) - sd(x)},
fun.max =function(x) {mean(x) + sd(x)},
geom = "errorbar")+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank()) +
facet_wrap(~class_3)
ggarrange(plotlist= list(angles_ms, angles_ms2, angles_ms3) , nrow= 3, ncol = 1, legend = "left")## Warning: Removed 5 rows containing missing values (geom_segment).
## Warning: Removed 5 rows containing missing values (geom_segment).
## Warning: Removed 5 rows containing missing values (geom_segment).
## Warning: Removed 5 rows containing missing values (geom_segment).
## Warning: Removed 5 rows containing missing values (geom_segment).
## Warning: Removed 5 rows containing missing values (geom_segment).
ggarrange(plotlist= list(degree_ms, degree_ms2, degree_ms3) , nrow= 3, ncol = 1, legend = "left")## Warning: Removed 1 rows containing missing values (geom_segment).## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_segment).
#Resumimos la data en el grafico tipo boxplot para evaluar los quartiles y su relación con valores médias:
ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure, fill = type_of_measure)) +
geom_boxplot() + stat_summary(fun = mean,
shape = 18,
size = 1) +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())## Warning: Removed 5 rows containing missing values (geom_segment).
ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure), color = "black") +
geom_boxplot(alpha = 0.3, aes(fill = type_of_measure)) +
stat_summary(fun.y = mean,
geom = "point",
shape = 18, size = 4,
aes(color = type_of_measure))+
facet_wrap(~class_2)+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())## Warning: `fun.y` is deprecated. Use `fun` instead.
ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure), color = "black") +
geom_boxplot(alpha = 0.3, aes(fill = type_of_measure)) +
stat_summary(fun.y = mean,
geom = "point", shape = 18,
size = 4,
aes(color = type_of_measure)) +
facet_wrap(~class_3)+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())## Warning: `fun.y` is deprecated. Use `fun` instead.
ortho_df2 %>% ggplot(aes(y= degree_spondylolisthesis,
x = 0)) +
geom_boxplot(fill = "red",
alpha = 0.4) +
stat_summary(fun = mean, shape = 18,
size =1)+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())## Warning: Removed 1 rows containing missing values (geom_segment).
ortho_df2 %>% ggplot(aes(x =0,
y= degree_spondylolisthesis,
fill= class_2),
color = black) +
geom_boxplot(alpha = 0.6,
aes(fill = class_2)) +
stat_summary(fun.y = mean,
geom = "point",
shape = 18, size = 5,
aes(color = class_2))+
facet_wrap(~class_2)+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())## Warning: `fun.y` is deprecated. Use `fun` instead.
ortho_df2 %>% ggplot(aes(x = 0, y= degree_spondylolisthesis), color = black) +
geom_boxplot(alpha = 0.4, aes(fill = class_3)) +
stat_summary(fun.y = mean, geom = "point", shape = 18, size = 5, aes(color = class_3))+ facet_wrap(~class_3)+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())## Warning: `fun.y` is deprecated. Use `fun` instead.
# Normality tests: Ahora vamos a comprobar si la distribución de las variables realmente sigue la distribución normal (hipotesis 0). Para evaluarlo vamos a utilizar el test de Kolmogorov-Smirnov, y la visualización con QQ-plot
# Kolmogorov-Smirnov
# Para realizar este test sobre el conjunto de data vamos a utilizar función map - una varaante de lapply() dentro del paquete tidyverse (libreria purr)
ortho_df %>%
select(where(is.numeric))%>%
colnames() %>%
set_names() %>%
map(~ lillie.test(ortho_df[,.])) %>%
map_dfr(., tidy, .id = "variable")## # A tibble: 6 × 4
## variable statistic p.value method
## <chr> <dbl> <dbl> <chr>
## 1 pelvic_incidence 0.0708 7.33e- 4 Lilliefors (Kolmogorov-Smirnov) n…
## 2 pelvic_tilt 0.0829 2.37e- 5 Lilliefors (Kolmogorov-Smirnov) n…
## 3 lumbar_lordosis_angle 0.0687 1.25e- 3 Lilliefors (Kolmogorov-Smirnov) n…
## 4 sacral_slope 0.0578 1.43e- 2 Lilliefors (Kolmogorov-Smirnov) n…
## 5 pelvic_radius 0.0593 1.06e- 2 Lilliefors (Kolmogorov-Smirnov) n…
## 6 degree_spondylolisthesis 0.173 5.13e-25 Lilliefors (Kolmogorov-Smirnov) n…# Vamos a mirar si el hecho de quitar el outlier (valor 418 del grado de espondilolistesis nos afecta al tipo de distribución)
ortho_df3 %>%
select(where(is.numeric))%>%
colnames() %>%
set_names() %>%
map(~ lillie.test(ortho_df[,.])) %>%
map_dfr(., tidy, .id = "variable") # Parece que no afecta, así que procedemos a seguir con el análisis con el data.frame completo.## # A tibble: 6 × 4
## variable statistic p.value method
## <chr> <dbl> <dbl> <chr>
## 1 pelvic_incidence 0.0708 7.33e- 4 Lilliefors (Kolmogorov-Smirnov) n…
## 2 pelvic_tilt 0.0829 2.37e- 5 Lilliefors (Kolmogorov-Smirnov) n…
## 3 lumbar_lordosis_angle 0.0687 1.25e- 3 Lilliefors (Kolmogorov-Smirnov) n…
## 4 sacral_slope 0.0578 1.43e- 2 Lilliefors (Kolmogorov-Smirnov) n…
## 5 pelvic_radius 0.0593 1.06e- 2 Lilliefors (Kolmogorov-Smirnov) n…
## 6 degree_spondylolisthesis 0.173 5.13e-25 Lilliefors (Kolmogorov-Smirnov) n…### QQplot
ggqqplot(ortho_df2, "angles", facet.by = "class_3", color = "type_of_measure")
ggqqplot(ortho_df2, "angles", facet.by = "class_2", color = "type_of_measure")
ggqqplot(ortho_df2, "degree_spondylolisthesis", color = "class_2", facet.by = "class_2")
ggqqplot(ortho_df2, "degree_spondylolisthesis", color = "class_3", facet.by = "class_3")
### Homogenidad - vamos a mirar si la data es homogénea
ortho_df %>%
select(where(is.numeric))%>%
colnames() %>%
set_names() %>%
map(~ fligner.test(ortho_df[,.] ~ class_2, data = ortho_df)) %>%
map_dfr(., tidy, .id = "variable")## # A tibble: 6 × 5
## variable statistic p.value parameter method
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 pelvic_incidence 18.9 1.40e- 5 1 Fligner-Killeen test of…
## 2 pelvic_tilt 14.0 1.80e- 4 1 Fligner-Killeen test of…
## 3 lumbar_lordosis_angle 25.3 4.98e- 7 1 Fligner-Killeen test of…
## 4 sacral_slope 16.4 5.11e- 5 1 Fligner-Killeen test of…
## 5 pelvic_radius 9.29 2.30e- 3 1 Fligner-Killeen test of…
## 6 degree_spondylolisthesis 89.8 2.58e-21 1 Fligner-Killeen test of…ortho_df %>%
select(where(is.numeric))%>%
colnames() %>%
set_names() %>%
map(~ fligner.test(ortho_df[,.] ~ class_3, data = ortho_df)) %>%
map_dfr(., tidy, .id = "variable")## # A tibble: 6 × 5
## variable statistic p.value parameter method
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 pelvic_incidence 7.21 2.72e- 2 2 Fligner-Killeen test of…
## 2 pelvic_tilt 27.3 1.19e- 6 2 Fligner-Killeen test of…
## 3 lumbar_lordosis_angle 19.3 6.56e- 5 2 Fligner-Killeen test of…
## 4 sacral_slope 7.28 2.63e- 2 2 Fligner-Killeen test of…
## 5 pelvic_radius 22.9 1.07e- 5 2 Fligner-Killeen test of…
## 6 degree_spondylolisthesis 125. 6.04e-28 2 Fligner-Killeen test of…Anova test:
Tras realizar varios tests de normalidad se confirma que la distribución de data NO es normal, asi que procedemos a realizar tests ANOVA, con la idéa en mente de que la distribución puede ser insuficientemente normal para conseguir resultados fiables. El test de homogenidad de data sugiere NO-homogenidad de data.
# Realizoms el test de anova para diferentes variables.
anova_PI <- aov(pelvic_incidence~class_3, ortho_df)
anova_S <- aov(degree_spondylolisthesis~class_3, ortho_df)
anova_PR <- aov(pelvic_radius~class_3, ortho_df)
anova_LLA <- aov(lumbar_lordosis_angle~class_3, ortho_df)
# La relación de las medias entre gurpos para la incidencia pélvica
par(mfrow=c(2,2))
plot(anova_PI)
par(mfrow=c(1,1))
plot(TukeyHSD(anova_PI))
# La relación de las medias entre gurpos para el grado de espondilolistesis
par(mfrow=c(2,2))
plot(anova_S)
par(mfrow=c(1,1))
plot(TukeyHSD(anova_S))
# La relación de las medias entre gurpos para angulo del radius de pelvis
par(mfrow=c(2,2))
plot(anova_PR)
par(mfrow=c(1,1))
plot(TukeyHSD(anova_PR))
# La relación de las medias entre gurpos para el angulo del lordosis lumbar
par(mfrow=c(2,2))
plot(anova_LLA)
par(mfrow=c(1,1))
plot(TukeyHSD(anova_LLA))
Los tests sugieren que la diferencia de los valores medias de algunas variables se difiere de forma significativa entre algunos grupos. (ver ejercicio 8)
#evaluamos un scatter plot para evaluar presencia de algun tipo de aglomeración.
ortho_df %>% ggplot(aes(lumbar_lordosis_angle, pelvic_incidence, color = class_3)) +
geom_point()
ortho_df %>% ggplot(aes(lumbar_lordosis_angle, pelvic_incidence, color = class_2)) +
geom_point()
# Miramos si se puede identificar algún numero de cluster
fviz_nbclust(ortho_df[,c(1,3)], kmeans, method = "wss")
fviz_nbclust(ortho_df[,c(1,3)], kmeans, method = "silhouette")
# Probamos classificar las valores segun la distancia entre los puntos del scatter plot para entre 2 y 3 clusteres
cluster2 <- kmeans(ortho_df[, c(1,3)], 2, nstart = 25)
table(cluster2$cluster, ortho_df$class_2)##
## Normal Abnormal
## 1 16 116
## 2 84 94summary(cluster2)## Length Class Mode
## cluster 310 -none- numeric
## centers 4 -none- numeric
## totss 1 -none- numeric
## withinss 2 -none- numeric
## tot.withinss 1 -none- numeric
## betweenss 1 -none- numeric
## size 2 -none- numeric
## iter 1 -none- numeric
## ifault 1 -none- numericfviz_cluster(cluster2, data = ortho_df[, c(1,3)])
cluster2 # Se evidencia que la classificación entre 2 clusteres es bastente complicada y la busqueda de algomeraciones por distancia no ha sido muy eficaz, y se acierta solo en 61.4% de casos.## K-means clustering with 2 clusters of sizes 132, 178
##
## Cluster means:
## pelvic_incidence lumbar_lordosis_angle
## 1 76.19685 68.75039
## 2 48.85381 39.45807
##
## Clustering vector:
## [1] 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [112] 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 2 2 2
## [149] 2 2 2 2 1 2 2 1 2 2 2 1 2 2 2 1 1 1 1 2 2 2 2 2 1 1 2 2 1 1 1 2 2 1 1 1 1
## [186] 2 1 2 1 2 1 1 2 1 2 1 1 2 2 2 2 1 2 1 1 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1
## [223] 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 43253.08 33224.06
## (between_SS / total_SS = 61.4 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"cluster3 <- kmeans(ortho_df[, c(1,3)], 3, nstart = 25)
table(cluster3$cluster, ortho_df$class_3)##
## Normal Hernia Spondylolisthesis
## 1 38 12 67
## 2 58 48 13
## 3 4 0 70summary(cluster3)## Length Class Mode
## cluster 310 -none- numeric
## centers 6 -none- numeric
## totss 1 -none- numeric
## withinss 3 -none- numeric
## tot.withinss 1 -none- numeric
## betweenss 1 -none- numeric
## size 3 -none- numeric
## iter 1 -none- numeric
## ifault 1 -none- numericfviz_cluster(cluster3, data = ortho_df[, c(1,3)])
cluster3 # y la clasterización en 3 clusteres se acierta en 74.1% de casos.## K-means clustering with 3 clusters of sizes 117, 119, 74
##
## Cluster means:
## pelvic_incidence lumbar_lordosis_angle
## 1 63.58681 53.46109
## 2 43.99372 35.12170
## 3 82.14937 76.54268
##
## Clustering vector:
## [1] 3 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [38] 2 2 2 2 2 2 2 2 2 2 1 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 1 1 2 2 2 2 2 1 1 2 1 2 2 1 2 2 2 2 2 2 2 2
## [112] 1 2 2 2 1 2 2 2 1 2 1 2 2 2 2 2 2 2 2 2 1 2 2 1 2 1 1 1 2 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 3 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 3 1 1 1 1 1 1 1 3 3 1 1 1 3 3
## [223] 1 1 1 1 3 1 3 3 1 3 1 1 3 1 1 1 3 3 1 3 1 1 3 1 1 1 3 3 3 3 1 3 3 1 3 3 3
## [260] 3 3 1 1 3 3 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 3 3 1
## [297] 3 3 1 3 3 3 3 3 3 3 3 3 3 3
##
## Within cluster sum of squares by cluster:
## [1] 17597.78 12236.25 21400.69
## (between_SS / total_SS = 74.1 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"# otra forma de realizar analisis de cluster en estos casos es clusterización jerárquica
cluster <- hclust(dist(ortho_df), method="complete")## Warning in dist(ortho_df): NAs introduced by coercionplot(cluster)
pltree(agnes(ortho_df, method = "complete"))
Ejercicio 8:
- A partir de los datos de origen y el estudio realizado (incluyendo todos los puntos anteriores), presentar un apartado de conclusiones. Esta sección debe incluir un resumen de los principales resultados obtenidos en apartados anteriores, que ayuden al lector a comprender el ámbito de estudio. Además, se valorará positivamente la coherencia de resultados y las justificaciones presentadas.
Mediante este estudio hemos podido resolver las dos preguntas principales hechas en el principio del trabajo.
Pregunta 1: ¿Existen diferencias significativas entre las variables angulares del grupo “Normal” respecto al grupo “Abnormal”?.
Según la muestra analizada, sí que existe una diferencia significativa entre algunas de las variables según los distintos grupos. Primeramente hemos visto cómo las medias de los pacientes con espondilolistesis para la incidencia pélvica y el ángulo de lordosis lumbar son superiores que en sujetos del grupo “Normal”. Este hecho denota que el desliz de las vértebras de los pacientes con espondilolistesis afecta a estas dos variables métricas más que a los otros sujetos. Por lo tanto, los resultados sugieren que la incidencia pélvica y el ángulo de lordosis lumbar podrían llegar a ser dos variables definitorias para la espondilolistesis. Para los pacientes con “Hernia”, no obstante, no había diferencias significativas respecto a los sujetos del grupo “Normal”. Por lo tanto, el desplazamiento vertebral (espondilolistesis) influye claramente más que la hernia en las variables métricas posturales.
Por otro lado, tenemos que no parecen existir diferencias significativas entre los grupos para el radio pélvico, la pendiente sacral y la inclinación pélvica, por lo que éstas podrían ser variables menos explicativas para las afectaciones estudiadas. De hecho, el grupo de pacientes con Hernia no difiere significativamente del grupo Normal para ninguna variable, dando a entender que las variables métricas estudiadas en este análisis no serían buenos distintivos para esta patología. Quizás con más estudios de muestras mayores y más representativas se podría contrastar los datos de las dos patologías y verificar si las variables utilizadas en estos ejercicios son buenas variables tanto definidoras como métricas.
Pregunta 2: ¿Existe alguna relación o correlación entre las variables del estudio?
Según los resultados de nuestros análisis, sí que existe una relación entre algunas de las variables.
Primeramente, hemos detectado una correlación significativa entre la incidencia pélvica con la pendiente sacral. Esta correlación se ha detectado mediante un análisis de regresión lineal y ha sido confirmado con un análisis de regresión múltiple. Con éste último también se ha detectado una correlación elevada entre la incidencia pélvica y la inclinación pélvica. La relación entre estas tres variables ha sido verificada con los datos en la bibliografía. Además, esta relación tiene sentido a nivel anatómico, ya que las modificaciones en la inclinación pélvica y/o la pendiente sacral hacen que el ángulo de incidencia pélvica cambia cuando se progresa la enfermedad (presencia de hernias o espondilolistesis).
Durante el análisis de los valores médias entre los grupos (ANOVA) se evidencia que, aunque la distribución de las variables no es normal, lo que hace que el test elegido no es muy fiable, en algunos grupos el valor média es diferente de forma significativa.
| Variable | Relacion |
|---|---|
| pelvic_incidence | NO varia entre el grupo normal y la hernia, pero varia en comparación con el grupo espondilolistesis |
| degree_spondylolisthesis | NO varia entre el grupo normal y la hernia, pero varia en comparación con el grupo espondilolistesis |
| pelvic_radius | NO varia entre el grupo espondilolistesis y la hernia, pero varia en comparación con el grupo espondilolistesis y el grupo normal |
| lumbar_lordosis_angle | Varia de forma importante entre los tres grupos. |
Durante el analisi de cluster se realiza 2 tipos de análisis - no jerárquico utilizando valores de distancia entre puntos, y jerárquica. Al realizar el test kmeans y al evaluar el data.frame se evidencia que el numero optimo de clusters se situa entre 2 y 3, por lo que se realiza analisis para las dos situaciones, objetivando mejor calidad de aciertamiento (74%) de clusterización para la situación con 3 clusteres.